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We present an extensive numerical study of the critical behavior of dimer models in three di- 
mensions, focusing on the phase transition between Coulomb and crystalline columnar phases. The 
case of attractive interactions between parallel dimers on a plaquette was shown to undergo a conti- 
nuous phase transition with critical exponents close to those of the 0{N) tricritical universality 
class, a situation which is not easily captured by conventional field theories. That the dimer model 
is exactly fine-tuned to a highly symmetric point is a non trivial statement which needs careful 
numerical investigation. In this paper, we perform an extensive Monte Carlo study of a generalized 
dimer model with plaquette and cubic interactions and determine its extended phase diagram. We 
find that when both interactions favor alignment of the dimers, the phase transition is first order, 
in almost all cases. On the opposite, when interactions compete, the transition becomes continuous, 
with a critical exponent rj ~ 0.2. The existence of a tricritical point between the two regimes is 
confirmed by simulations on very large size systems and a flowgram method. In addition, we find 
a highly-degenerate crystalline phase at very low temperature in the frustrated regime which is 
separated from the columnar phase by a first order transition. 



PACS numbers: 05.30.-d, 02.70.Ss, 64.60.-i, 75.10.-b 



I. INTRODUCTION 

Originally suggested^ as descriptive of adsorption of 
molecules on a substrate (a motivation that has been 
renewed in recent experiments^), dimer models have at- 
tracted the interest of researchers in various branches of 
physics, ranging from statistical and condensed matter 
physics to high-energy physics'^. Their distinctive pro- 
perties essentially result from the close-packing condi- 
tion, which imposes that on a lattice, each site should be 
part of one and only one dimer. This strict condition ge- 
nerates strong correlations between degrees of freedom, 
even when interactions are absent from the system. 

Classical dimer models have been originally studied in 
statistical mechanics, with the famous result that non- 
interacting dimer models on planar graphs can be sol- 
ved exactly using PfafRans'*'^. On the square lattice for 
instance, it was shown using subsequent techniques that 
dimer correlation functions decay algebraically with dis- 
tance^. It was latter shown that dimer models can also 
be viewed as dual versions of Ising models^ and generate 
the ground-state manifolds of fully frustrated Ising mo- 
dels''. For bipartite lattices in three dimensions, dimers 
can be represented by an effective magnetic field living on 
the bonds of the lattice. The close-packing condition for 
the dimers encodes a Gauss law for the magnetic field*^. 
Postulating a quadratic dependence on this field of the 
effective cntropic action suggests the existence of dipo- 
lar dimcr-dimer correlations. Monte Carlo simulations^ 
indeed confirm this picture with a great accuracy, and 
the corresponding phase of dimers on 3d bipartite lat- 



tice is often referred to as a Coulomb phasc^ with this 
electromagnetic analogy in mind. 

In some sense, dimer models form the "Ising model" 
of local constraint, spelling out their ubiquity in physics. 
In condensed-matter, dimer models emerged as classical 
counterparts of quantum dimer models^" (which proper- 
ties for particular values of parameters are determined by 
those of the classical problem), as well as effective mo- 
dels for magnetization plateaus in frustrated magnets^ ^. 
Dimer models show also close similarities with spin-ice 
systems^^ which can also host a Coulomb phase. There 
the close-packing condition translates into the "ice rule" . 

Intriguing physics take place in classical dimer models 
when interactions are present. The perhaps simplest case 
to study is to add local interactions which favor paral- 
lel alignment of dimers on a plaquette of the lattice. In 
two dimensions on bipartite lattices, the system under- 
goes a phase transition from a columnar phase at low 
temperature to a disordered critical phase at high tem- 
perature^'^. The transition is of the Kosterlitz-Thouless 
type. It is possible to obtain a field theoretical descrip- 
tion of the transition in terms of an height model which 
predicts accurately the behavior of the correlation func- 
tions of different observables^^"^^. The situation is much 
less clear in three dimensions. At high temperature on 
the cubic lattice, the system is located in the Coulomb 
phase which is destabilized as the temperature is lowered 
towards a columnar order of dimers. Quite surprisingly, 
the transition between the critical Coulomb phase and 
the ordered phase is continuous^^ . Critical exponents es- 
timated from the numerical simulations do not appear to 
be those of a known "simple" universality class, although 



being very close to those of a tricritical theory. 

A field theoretical description of the transition obser- 
ved in the 3d classical dimer model is not easy. In parti- 
cular, it cannot be properly addressed in the traditional 
Ginzburg-Landau formalism since the correlations bet- 
ween degrees of freedom in the disordered phase decay 
algebraically and not exponentially. Different attempts 
to provide for a field theoretical description of the criti- 
cal point have led to a representation in the continuum 
in terms of two complex matter fields coupled to a non 
compact C/(l) gauge field^^"^^, a theory known as the non 
compact CP^ (NCCP^) theory. The problem is that it 
is not clear at present if this theory possesses in fact an 
infra-red fixed point. Efforts to simulate lattice versions 
of the NCCP^ model either lead to a weakly first-order 
transition^*^ , or to a continuous transition with uncon- 
ventional critical exponents^^'^^. Of course, it is possible 
that the NCCP^ theory possesses a tricritical point and 
that the different simulated microscopic models all cor- 
respond to the same theory but fiow in the continuum 
toward different parameter regimes. 

Another possibility is that the microscopic dimer 
model of Ref. 16 sits at a tricritical point, which would 
be the essentially unique way of reconciling the observed 
continuous transition with a Ginzburg-Landau approach 
based on symmetry-breaking. Note that the value of the 
critical exponents a ~ 0.5, v ~ 0.5 and rj ~ 0^^ are 
suggestive of this scenario as they are close to the ones 
of a 0{N) tricritical theory. This putative tricritical 
point could be the one of the NCCP^ theory, or of 
another field theory yet to be specified. It is however 
unclear why the dimer model should be precisely 
located at a tricritical point : this usually requires a 
fine-tuning of parameters, and there is no other parame- 
ter than temperature in the original microscopic model'^^. 

In this paper, we step aside from field theoretical 
considerations and instead provide new valuable data 
regarding the critical behavior of 3d classical dimer 
models. We have carried an extensive Monte Carlo (MC) 
simulation of a dimer model consisting of a four-dimers 
interaction on a cubic lattice, in addition to the usual 
attractive two-dimers plaquette interactions. The cubic 
interaction corresponds to a coupling between four 
parallel dimers sitting on the edges of a cube which 
can be attractive (non frustrated regime) or repulsive 
(frustrated regime). The system is studied with a worm 
MC algorithm which allows us to sample systems of 
linear size up to L = 180. Our results suggest that the 
phase transition between Coulomb and columnar phases 
is first order almost everywhere in the non frustrated 
regime, and second order in the frustrated side - forcing 
the existence of a tricritical point in between these two 
regimes. Moreover, we find that the critical exponents in 
the strongly frustrated regime are different from the ones 
measured in the absence of the cubic interaction, with 
rj ~ 0.2 and v '^ 0.6. In order to rule out a very weak first 
order transition for all the range of parameters (which 



is always possible), we performed a flowgram analysis 
which clearly indicates two collapses above and below 
the tricritical point. Finally, at very low temperature 
in the strongly frustrated regime, we observe a new 
crystalline phase resulting from the destabilization of the 
columnar phase. This phase has a degeneracy growing 
extensively with the linear system size and is separated 
from the columnar phase by a first order transition. 
The final phase diagram that we obtain is presented 
in Fig. 1. We note that similar results were recently 
obtained by Papanikolaou and Betouras^^ in a different 
deformation of the dimer model. A comparison between 
the two works will be made in Sec. VI. 



The plan of the paper is the following. In Sec. II, we 
describe the model, the algorithm and the relevant ob- 
servables for its study. In Sec. Ill, we present our results 
in the non frustrated regime, when both interactions are 
attractive. In Sec. IV, we analyze the frustrated regime 
where the plaquette interaction is attractive but the cu- 
bic interaction is repulsive. The final form of the phase 
diagram is obtained with a flowgram analysis, which we 
present in Sec. V. We finally discuss the implications of 
our results in Sec. VI. 
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Figure 1: Phase diagram of the attractive dimer model with 
cubic interactions. Dots indicate values of interactions where 
simulations were performed. Dashed lines and green dots re- 
present second order transitions. Solid lines and red dotes cor- 
respond to first order transitions. The white dot represents a 
point where the nature of the transition is still unclear. 



II. DEFINITIONS 



A. System 



The system we consider is a cube of linear dimension 
L (total number of sites N = L^) covered by hard-core 
dimers. Only dimer configurations C obeying the close- 
packing condition contribute to the partition function. 
The partition function reads : 



Z = Y, cxp{~/3Ec 



(1) 



and the energy of an allowed configuration C is given by 



Ec - V2{Nii +N^+Nfj)+ W4iVeubcs, 



(2) 



where the first term is proportional to the number of 
plaquettes in the configuration containing two parallel di- 
mers (referred in the following as "parallel plaquettes") 
and the second counts the number of unit cubes sustai- 
ning four parallel dimers (see Fig. 2). Both occurences 
of terms in a given dimer configuration are illustrated in 
Fig. 3. In the remainder of this study, we will restrict 
ourselves to the case of attractive plaquette interactions 
V2 < while cubic interaction V4 can be attractive or re- 
pulsive. We investigate the properties of the system as a 
function of the ratio x — W4/W2 and temperature T = 1/(3. 
When a; > 0, both interactions have the same sign and 
the system will be said to be non frustrated, as both 
interactions favor the same columnar ground-states. On 
the opposite, when a; < 0, the system will be in a frus- 
trated regime (v2 and V4 terms compete). In general, we 
take W2 = — 1 and vary V4. The only exception is the two 
limits X — > ±00 where we consider the system in absence 
of plaquette interactions {v2 = 0^) and take W4 = ±1. 




Figure 2: The three different possibilities with four parallel 
dimers parallel on a unit cube. Each pattern contributes V4 
to the energy Ec in Eq. 2. 



We simulate the dimer model by means of a wornr 
Monte Carlo algorithm with a local heat-bath detailed 
balance condition^^. Compared to a local Metropolis al- 
gorithm, autocorrelation times are drastically reduced 
with a worm algorithm due to the use of non-local moves, 
which allows us to reach systems of linear size up to 




Figure 3: Illustration of occurences of the V2 (top) and V4 
(bottom) terms of the model defined in Eq. 2 : interacting 
plaquettes (top) and cube (bottom) are represented with sha- 
ded surfaces. 



L = 180. As we will see in the next section, the abi- 
lity to simulate very large system sizes is of crucial im- 
portance to distinguish between continuous and weakly 
first-order transitions. We would indeed like to empha- 
size that for systems like dimer models, or others which 
contain a priori long-range correlations, one should be 
particularly cautious regarding the issue of finite-size sca- 
ling. Here, for the largest system sizes, up to 5.10^ sweeps 
have been carried (we define one sweep by performing en- 
ough worm moves such that on average every site of the 
lattice is visited). The convergence of simulations is che- 
cked by looking at autocorrelation times of the various 
observables, as obtained from a binning analysis. 



B. Observables 

The observables in our system are of three kinds : a first 
group is made of the thermodynamic quantities such as 
the average energy, a second type of observables is related 
to the columnar ordering, and the third type is related to 
the stiffness of the system (fluctuations of dimcr fluxes). 



1. Thermodynamic quantities 

A phase transition can generally be detected by mo- 
nitoring the probability distribution of the energy of the 
system. For a first order transition, the average energy 
{E) is discontinuous and exhibits a latent heat when the 
temperature is lowered. For a second order transition, the 
average energy is continuous but its first derivative, the 
specific heat per site C^/N, obeys the scaling law : 



~N N dT 



1 d{E) _ (E^) - {Ef 



T^L^ 



Cl''^ + AL°'/\ (3) 



The first term C^°^, corresponding to the regular part of 
the specific heat at criticality, is often forgotten in fits 
of numerical data as the divergence of the second term 
usually dominates (for a > 0). However, at several points 
of the phase diagram, we find that this term cannot be 
neglected as the divergence of the specific heat can turn 
quite slow. In this situation, one can take C™^ as either 
a new fitting parameter or as given by results obtained 
on small systems where the second term is negligible. 
In these cases, we take a conservative approach for the 
error bar on a/v and quote a result which encloses all 
possibilities. 

In comparison, for a first order transition, the specific 
heat diverges like the volume : Cy/N ex L^. Finally, ano- 
ther mean to distinguish between first and second order 
transitions is to consider the whole histogram of energy 
at the transition temperature, as obtained in the Monte 
Carlo simulation. We expect for a first order phase tran- 
sition the appearance of double peaks, centered at the 
average energy of the two co-existing phases. These peaks 
will appear only for samples with size above (or close 
to) the correlation lenght at the transition, and should 
keep being separated when increasing system size. For 
a second order phase transition, the histogram should 
contain a unique peak. 



2. Columnar order parameter 

The local columnar order parameter m(r) is defined 
with respect to the dimcr occupation number at each 
site n(r) : 



m(r) = (-l)'"n(r). 



(4) 



The global order parameter reads C = ^ J^r ^^i^) '"^^^ i^ 
normalized such that the six columnar states correspond 



to : C = {±1, 0, 0}, {0, ±1, 0}, {0, 0, ±1} and C = (|C|) = 
1. One also considers the corresponding susceptibility x- 
In particular, for a second order transition : 



X/N = 



jC) (ici: 

i3 



L^-\ 



(5) 



while for a first order transition x/-^ ^ ^^- Finally, the 
Binder cumulant : 



B = 



(6) 



is a scale-invariant quantity in the case of a continuous 
transition, and should thus exhibit a crossing point at 
criticality as a function of the system size. Moreover, the 
finite size-scaling (FSS) of its derivative with respect to 
the temperature : 



dB 
dT 



ex L 



Iju 



(7) 



gives a direct access to the critical exponent v. 



3. Stiffness 

The (inverse) stiffness encodes fluctuations of dimcr 
fluxes across a plane^'^^ : 



^-'- E 



a—x,y,z 



ZL 



(8) 



where the flux (pa is the algebraic number of dimers cros- 
sing a plane perpendicular to the unit vector a . Given a 
lattice direction, the contribution to the flux is +1 for a 
dimer going from one sublattice to the other and —1 for 
the reverse situation. The stiffness is finite in the Cou- 
lomb phase, refiecting the presence of dipolar correlations 
between dimers. On the other hand, the columnar phase 
is robust to insertion of fluxes and K'^ vanishes expo- 
nentially with system size. At a second order phase tran- 
sition, the quantity LK^^ should be scale invariant^^ and 
the scaling of its derivative : 



dK- 



dT 



^l/u 



(9) 



provides another access from the high temperature side 
to the exponent v. In general, the error bars that we 
quote on exponents include at the same time errors due to 
the fitting procedure (which we measure by considering 
stability of fits with exclusion of a few data points) , errors 
from the determination of critical temperature as well as 
statistical errors. 



III. NON FRUSTRATED SIDE : V2<Q,Vi< 0. 

We start the discussions of our numerical results on 
the non-frustrated side W4 < and f 2 < of the phase 



diagram. When both interactions are attractive, we ex- 
pect to find the same phases as in the simple attractive 
plaquette model : a six-fold degenerate columnar phase 
at low temperature and a Coulomb phase with dipolar 
correlations at high temperature. 

At first, let us consider the extreme case where only 
attractive cubic interactions are present {v2 = or 
tanh(w4/w2) = 1 in the phase diagram of Fig. 1). The 
evolution of the energy, columnar order parameter C and 
inverse stiffness K~^ for small system sizes is presented 
on Figure 4. As expected, the columnar order is non- 
zero only at low temperature where the inverse stiffness 
vanishes. But in constrast with the attractive plaquette 
model, both quantities exhibit a strong discontinuity at 
the transition temperature Tc ^ 1.1. In fact, the energy 
also displays such a jump, characteristic of a latent heat. 
This shows undoubtedly a first order transition. 
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Figure 4: Left : Evolution of the columnar order parameter 
C and inverse stiffness K~^ as a function of temperature for 
V2 = and V4, = —1. Right : Temperature dependence of the 
average energy per site (E)/N. Here L = 16. 



To settle definitively the nature of the transition as well 
as to benchmark the method, we also study the histogram 
of energy observed during the simulation, which corres- 
ponds to the probability distribution of the energy P{E). 
The appearance of a double peak distribution at the criti- 
cal temperature is a typical sign of a first order transition. 
For V4^/v2 = oo, we can easily detect this double peak for 
system sizes as small as L = 8 (see Fig. 5). 

We now introduce a small attractive plaquette interac- 
tion V2 and repeat the procedure by tracking the tem- 
perature at which columnar order sets in and inverse 
stiffness vanishes. We observe that the nature of the 
transition remains discontinuous but that the correla- 
tion length grows as the ratio U4/W2 is decreased : for 
Vi/v2 = 0.8 we find double peaks only at sizes L > 16, 
for i'4/w2 = 0.6 at sizes L > 32 and for ■i;4/i;2 = 0.4 at 
L > 80 (see Fig. 5). Finally, for systems close to the pure 
plaquette model W4/W2 — 0, it becomes extremely difficult 
to distinguish double peaks in the energy distribution. In 
fact, the histogram for W4/W2 = 0.2 shows a slightly de- 
formed single peak for L = 140. We were not able to see 
any sign of double peaks for vA/v2 = 0.1 up to L = 140. 
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Figure 5: Energy per site probability -P(e) as a function of 
energy per site e = E/N, for different ratios U4/W2 > at 
the critical temperature separating columnar and Coulomb 
phases. The size L corresponds to the minimal length for 
which the histograms start to display a double peak distri- 
bution. For v^jvi = 0.2, a deformed single peak is observed 
in the distribution for L = 140, the maximal sample size si- 
mulated for this parameter. 



Another possibility to discern a discontinuous transi- 
tion is to measure the critical scaling of the maximum of 
the susceptibility and the specific heat per site. At a first 
order transition, both quantities should diverge like the 
volume L^ . In terms of the critical exponents introduced 
in section II B, this would correspond to effective critical 
exponents ^| „ = 3 and 77off = —1. We have determi- 
ned the scaling laws of these two quantities for several 
values of the ratio W4/W2 (see table I) and find that when 
vaIv2 is large, the exponents agrees with the first order 
values. As we approach W4/U2 — 0, the scalings of Cy and 
X get smoother and the exponents closer to the values 
of Ref. 16 ^1 „ ~ 1 and 77eff ~ 0. At this point, it is 
not possible to conclude on the order of the transition at 
Vi/v2 — 0.1. The transition can be either continuous or 
very weakly first order. 
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Table I: Critical temperature and effective critical exponents 
for different non frustrating coupling ratios. For each ratio, the 
exponents have been measured by a FSS analysis using Lmin 
as the minimal size. The values at Vi/v2 = are taken from 
Ref. 16. 



IV. FRUSTRATED SIDE : ua < U4 > 

We now turn to the analysis of the system where pla- 
quette interactions are attractive but cubic interactions 
are repulsive. As the two interactions compete with each 
other, we can expect, at least at low temperature and in 
the regime W4/W2 ^ ^1, that new phases may appear. 

Let us again first discuss the extremal case where only 
cubic interactions are present {v2 = 0). We observe that 
the columnar order parameter vanishes for all tempera- 
ture while the inverse stiffness K~^ remains always non- 
zero and finite (see Fig. 6 left). Moreover, although the 
specific heat displays a maximum (see Fig. 6 right), it 
does not display any dependence on the size of the sys- 
tem and cannot be related to a critical phenomena. Thus, 
we obtain the surprising result that the system is always 
disordered when t;4 > and v^ =0. In other words, the 
Coulomb phase can accomodate for having no cubes oc- 
cupied by four dimers (as in Fig. 2), even down to very 
low temperatures. Even more, the Coulomb phase ap- 
pears to be strenghten in this situation as the inverse 
stiffness increases slightly as temperature is lowered. 



two critical temperatures get closer such that it becomes 
more difficult to detect the peak at Tc2- 
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Figure 6: For U4 = 1 and V2 — 0. Left : Evolution of the stif- 
ness K~^ and columnar order C as a function of temperature 
for L = 16. The non-zero value reached by C is due to the 
finite size of the sample. Right : The height of the peak of 
the specific heat per site does not display any dependence on 
system size. 



As soon as a finite attractive coupling V2 is introduced, 
we find that the specific heat per site displays two peaks 
as a function of temperature : one at a lower temperature 
Tcj^ which is strongly diverging with the system size and 
another at an upper temperature T^^ which diverges very 
slowly. Fig 7 and its insets display results at V4/V2 = — f , 
which are typical of what we observe in the frustrated re- 
gime. The first peak of the specific heat is associated with 
the freezing of the columnar order parameter to a value 
smaller than 1 at T < Tci , the second to the appearance 
of the Coulomb phase for T > Tc^ (see upper inset of fi- 
gure 7). In the next two sections, we will detail the nature 
of these two phase transitions and of the phase separa- 
ting them. We note that as the ratio W4/W2 — > —00, the 
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Figure 7: Evolution of the specific heat per site as a function 
of temperature for W4/i'2 = —1. The first peak at T^ ~ 0.5 is 
associated with the crystallisation of the system in an ordered 
phase different from the columnar phase, which can be readily 
seen by the evolution of the columnar order parameter (see 
upper inset). The second peak, at T^j ~ 0.95, delimitates a 
region with non-zero columnar order from a region with non- 
zero inverse stifl'ness. Upper inset : evolution of the inverse 
stifl'ness K~^ and of the columnar order parameter C for L = 
16. Lower inset : zoom of the specific heat close to Tc2 ■ 



A. Low temperature phase transition Tci 

We now discuss the nature of the phase below the lo- 
wer transition by first considering the evolution of the 
T = ground-states as a function of the frustration ra- 
tio V4/v2. For large positive values of V4 (but still with 
V2 — — f), we expect the columnar ground-states to be- 
come unstable as they maximize the number of parallel 
cubes. In this limit, we must find dimer configurations 
which have exactly zero parallel cubes but that can none- 
theless support a maximal number of parallel plaquettes. 
We find that there are several configurations (in fact an 
exponential number) that satisfy this frustrating condi- 
tion. For instance, we present in Fig. 8 on the example 
of a L = 6 cube three configurations satisfying these two 
constraints. Configuration A is made of a unit pattern 
consisting of two planes repeating each other. The unit 
pattern of configurations B and B' posseses three planes. 
Configurations B and B' are simply related by a unit 
translation of the two bottom planes. If we now consi- 
der larger systems, it is easy to see that we can use the 
same plane patterns at will, by randomly choosing B or 
B' every three planes and this, without creating parallel 
cubes. Therefore the degeneracy of the ground-state is 
at least growing like 2^'^, that is exponentially with the 
linear system size. It is possible that by forming other 
cost-free defects within a plane, the degeneracy is even 



higher ex e° , but we have not made any further inves- 
tigations in this direction. 
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Figure 8: Three configurations on the L — 6 cube which 
maximize the number of parallel plaquettes without having 
any parallel cubes. Configuration A is built out of a unit pat- 
tern containing two planes which repeats itself three times. 
This unit pattern possesses 42 parallel plaquettes, leading to 
a number of 7/12 parallel plaquettes per unit site. Confi- 
gurations B and B' both consist of a unit pattern of three 
planes repeating itself twice. The unit pattern possesses 63 
parallel plaquettes, leading also to a number of 7/12 parallel 
plaquettes per unit site. Configuration B' is obtained from 
configuration B by shifting the dimer pattern on the two lo- 
west planes by one unit cell on the right direction. Mixing B 
and B' patterns, one can create several dimer configurations 
having the same energy. 



For wliat thermodynamics is concerned, it is easy to 
check that all such configurations with no parallel cubes 
have on average 7/36 plaquettes that are parallel. We re- 
fer to the corresponding phase as Crystal II. While it may 
be possible to define correctly an order parameter for this 
crystal (and this in spite of the high ground-state dege- 
neracy), we simply concentrate on a simple comparison 
between the energy of the Crystal II ground-states with 
those of the columnar phase. A columnar ground-state 
satisfies 1/3 of the possible parallel plaquettes, and one 
parallel cube out of two. Since the number of plaquettes 
(cubes) is three times (equal to) the number of lattice 
sites, we obtain : 
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(10) 



That is, the Crystal II phase should be favored as soon as 
va/v2 < —10/12 ~ 0.833. We now compare this simple es- 



timate with numerical simulations. Considering the evo- 
lution of the columnar order parameter as a function of 
T for different ratios W4/f2, we find that C converges to 
1 at very low temperature for v^jv^ > —0.8, while for 
W4/W2 < —0.8 it converges towards a smaller value (see 
Fig. 9 top). This indicates that the lowest energy confi- 
gurations are no longer the columnar ones. A further in- 
dication of the Crystal II phase is given by the average 
number per site of parallel cubes {Ncuhc}/L^ and paral- 
lel plaquettes (A'^piaquctto)/^'^- As expected, we find (see 
Fig. 9 bottom) that {Ncnhc)/L^ decreases from 1/2 to 
and {Npii^quette) / L^ from 1 to 7/12 quite abruptly as soon 
as f4/'y2 ^ —0.8. Note that the estimate W4/u2|c — —0.8 
that we obtain is quite rough as it is affected by the 
chosen grid in V4/V2, the moderate size of the sample 
(L = 32) and the finite temperature used in our simula- 
tions. Given this, it can be considered as in good agree- 
ment with the exact value —10/12. 

The abrupt behaviour observed in Fig. 9 tends to in- 
dicate a first order transition between the Crystal II and 
columnar phases. This is confirmed by the strong diver- 
gence of the specific heat (Fig. 7), but also by the ap- 
parition of a latent heat (see Fig. 9 top). We find that 
in the full phase diagram, the transition at Tc^ is always 
of first-order nature. When i'4/w2 ?J —0.8, the Crystal II 
phase and the low-T phase transition disappear. 



B. High temperature phase transition Tc2 

The high temperature phase transition corresponds to 
the simultaneous emergence of dipolar correlations at 
high temperature and disappearance of the columnar or- 
der. Before performing a detailed scaling analysis, we al- 
ready make an important statement : in all simulations 
for V4 > 0, we found no evidence for a first-order phase 
transition at T^^ between the columnar and the Coulomb 
phases. This has been checked in all observables at hand 
(thermodynamics, related to columnar order or the Cou- 
lomb phase), including energy histograms at the transi- 
tion. We will come back to this issue at the end of this 
section, and in Sec. V. 



1. Strongly frustrated regime 

In order not to be influenced by any crossover effect, we 
first concentrate on the transition at Tc2 far away from 
the putative tricritical point at V4 = 0. For V4/V2 = —10, 
we performed large scale simulations and applied the FSS 
analysis to calculate the critical exponents. Our set of 
data is presented in Fig. 10. Let us first concentrate on 
the specific heat per site (see Fig. 10 top). At the transi- 
tion, it exhibits a very slow growth with system size, sug- 
gestive of a continuous transition. In fact, the divergence 
is so small that the regular part Cl°^ of the specific heat 
contributes the most even for L = 140. The best fit we 
obtained for system sizes ranging from L = 32 to L = 140 
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Figure 9: Top : Evolution as a function of temperature of 
the columnar order parameter C for different ratios Vi/v2 at 
L — 32. For V4/v2 < —0.8, C does not converge towards 1 
at T = 0. Inset : Evolution of the average energy per site. 
A discontinuity corresponding to the phase transition can ea- 
sily be deteted. Bottom : Evolution of C, {Npiaquette)/N and 
{Ncuhe)/N as a function of U4/«2 for T = 0.3 and L = 32. 



gives an estimate of a/v ^ 0.4 (Fig. 10 top inset). Un- 
fortunately, this estimate varies a lot when considering 
another subset of system sizes. For instance, discarding 
the point at L = 140 leads to an estimate a/iy — 0.6 
and discarding the two points L — 120 and L — 140 
leads to a/i^ ^ 0.8. Thus, while we cannot conclude on 
the precise value of a at this point, it seems at least to 
be quite small. The evolution of the Binder cumulant 
and of the product K^^.L is presented in the middle 
panel of Fig. 10. For both quantities, we observe a cros- 
sing point, in agreement with a second order transition. 
The two crossing points temperatures are very close : 
^Binder = 0.672 and Tstiff„css = 0.6715. The thermody- 
namic measurements of the derivative quantities dB/dT 
and LdK~^ /dT, shown in the insets, allows us to have ac- 
cess to two independent estimates of the exponent v. We 
find ^'Binder ~ ^'Stiffness = 0.63(4). Thcse cxpoucnts are 
compatible with the Ising and XY universality classes in 
Sd. Moreover, this is consistent with a small value of a 
assuming hyperscaling a — 2 — dv. The value of the Bin- 
der cumulant and the stiffness at the critical point. Be 
and {K~^ ■L)c, also furnish valuable information, because 
these quantities are also universal. Checking the crossings 
of the three largest system sizes of our system, we find 



1.28 < Bc< 1.30 and 0.16 < {K''^ ■ L)^ < 0.20. While 
the value of B^ is consistent with the result obtained 
for the pure plaquette model Bc{v4^/v2 = 0) = 1.27(1), 
the value at criticality for the stiffness is rather smaller 
{K-^ ■ L)c{vi/v2 = 0) = 0.28(2). Finally, we discuss the 
scaling of the columnar susceptibility (see Fig. 10 bot- 
tom), obtaining r] = 0.25(3). The fact that t] is large and 
positive for the strongly frustrated system is an unam- 
biguous result of our study, as it is robust for instance 
on the set of sizes chosen to perform the FSS analysis. 
Such a strong 77 rules out the possibility of a simple tran- 
sition such as the Ising or XY universality class. Finally, 
we note that all the exponents differ considerably from 
those measured in Ref. 16, implying a different type of 
transition. This can be already be seen at the qualitative 
level as the specific heat does not display any (strong) 
divergence. 



2. Medium and weakly frustrated regime 

We repeated the same analysis for different values of 
W4/U2 on the frustrated side, measuring the exponents a, 
1/ and 77. In particular, we carried out extensive simu- 
lations at V4^/v2 = —1 and W4/U2 = —0.2. Results are 
summarized in Tab. II. For vn/v2 = —1, the exponent 
V is compatible with the one obtained for W4/f2 — —10, 
and while estimations of 77 are slightly different, the ano- 
malous dimension is clearly non-zero in both cases. The 
ratio of exponents a/v is again the hardest to determine 
reliably due to the important contribution of the regular 
part in the specific heat. For Vi/v2 = — 1, C^/N has ho- 
wever a clear diverging tendency and it is then easier to 
extract a/v (see figure. 11). We find a/v ~ 0.35, a value 
in accordance with the rough estimate at vA/v2 = — 10 
and with the hyper-scaling relation : a — 2 — dv. In any 
case, the set of exponents we obtain for Vi/v2 = —10 and 
Vi/v2 = — 1 are clearly different from the one obtained 
for the model with no cubic term. For v^lv2 = —0.2, 
the values are on the contrary compatible with the ones 
obtained from Ref. 16 : a ~ 0.5 and 77 ~ 0. Interestin- 
gly, we find two non overlapping estimates of v from the 
Binder cumulant and the stiffness derivatives. That may 
indicate a possible crossover. As a further information, 
we also give the values of Binder cumulant and stiffness 
at criticality. 

The results on the frustrated regime bring some in- 
terrogations. Our results suggest that there are two dif- 
ferent sets of critical exponents (and therefore universa- 
lity classes) for the Coulomb-columnar phase transition 
in the extended dimcr model : one for the highly frustra- 
ted regime with a ^ 0.2, v ~ 0.63 and 77 ^ 0.2 and one 
close to the point v^ — Q where the exponents are close 
to those of the 0{N) tricriticality class. Moreover, in the 
non frustrated regime, we have seen in section III that 
the transition between Coulomb and columnar phases is 
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Table II: Critical exponents for different frustrating coupling ratios. For the ratio a/v, the exponents have been measured by a 
FSS analysis using sizes comprised between Lmin = 32 and Lmax = 140. For the other exponents, we can limit ourselves to larger 
system sizes between Lmin = 80 and Lmax = 140. In the last row, we remind the results for the pure plaquette model (taken 
from Ref. 16). The symbol * denotes one case where the determination of a/v is impossible due to the strong contribution of 
the regular part of the specific heat (see Sec. II). 



clearly first order, at least when v/^jvi > 0.4. A natural 
interpretation of these data is that there is a tricritical 
point for a value w| in the vicinity of W4 = 0, separating 
a continuous transition in the frustrated case to a dis- 
continuous one in the non-frustrated case. The presence 
of a tricritical point can influence the effective critical 
exponents measured for values of v^jvi in the vicinity 
of v\lv2- Despite the large samples used in our simula- 
tions, we are not able to precise the exact location of this 
tricritical point. In particular, there is no formal reason 
to believe that the tricritical point is located exactly at 
114 — 0. In the next section, we will show that the tricriti- 
cal point should at least be located on the non frustrated 
side, that is W4/W2 ^ 0. 

Another possible explanation that we cannot exclude 
o 'priori is that the transition is always first order. In 
fact, by looking at the evolution of the different expo- 
nents from the non frustrated to the frustrated side, 
one can perfectly imagine that the correlation length 
grows continuously but remains finite at any va^Jvi. That 
would mean in particular, that the exponents found in 
Ref. 16 for the pure plaquette model are artifacts caused 
by a very weakly first order transition. There are seve- 
ral examples of models where reports of unconventional 
continuous phase transitions have been made and which 
were finally found to be weakly first order. In order to 
rule out this scenario, we present in the next section a 
flowgram analysis combined with a study of histograms 
at very large system sizes. 



V. FLOWGRAM ANALYSIS 

The fiowgram method is an advanced version of the 
FSS analysis which was proposed by Kuklov and coau- 
thors^'*. Consider two points located on the critical line 
separating the Coulomb and columnar phases in Fig. 1. 
The method relies on the demonstration that the large 
scale behavior for one given point is identical to that of 
the second point where the nature of the transition can 
be easily determined. If this is true, then the two points 
are in the same critical regime and the nature of the tran- 
sition remains unchanged all between the two points. On 
the contrary, a change in the nature of the transition 
must occur if this is not true. 



The key elements of the method are to (i) introduce 
a definition of the critical point for finite-size systems 
consistent with the thermodynamic limit and which is in- 
sensitive to the transition order and (ii) compute a quan- 
tity Q which is scale invariant at criticality, vanishes in 
one phase and diverges in the other. To define the opera- 
tional critical temperature, we assume a finite probability 
of having a non-zero flux at criticality^*^ : 



1 - P(0 = 0) 



A, 



(11) 



where A is some constant which exact value is not rele- 
vant for the rest of the method. In practice, we chose A in 
such a way that it is close to typical values found at the 
transition for a large system size for W4/W2 = —0.2. We 
then consider Q = L.K~^, as this product is indeed scale 
invariant for a second order transition, vanishes in the 
columnar phase (as we expect K~^ to vanish exponen- 
tally with system size) there and diverges in the Coulomb 
phase (as K^^ is a constant). 

The flows Qvi/v2{Li) for several values of W4/f2 in the 
interval [—0.6, 0.7] are presented on Fig. 12 top. The flows 
can be roughly divided into two groups : for W4/W2 < 0.2, 
the flows show a very slow divergence with the system 
size. For W4/U2 > 0.2, the flows are strongly diverging. 
The second group of flows is thus associated with the 
strongly discontinuous transition. At this stage we can- 
not draw any conclusion about the first group of curves 
since it might be that all curves diverge and are actually 
connected by a scaling transformation. By this, we mean 
that we need to check if there exists no renormalization 
fonction g{v4^/v2) such that plotted as a function of the 
renormalized length : 



Loff = g{vi/v2)L, 



(12) 



all the different flows collapse into a single master curve : 

Qvi/v2=a{Lcs) ^ Qvi/v2=b{Lcfi) y ia,b). (13) 

To search for such a transformation, we first try to col- 
lapse the flows two by two, starting from the largest po- 
sitive values of Vi/v2- Fixing g(vi/v2 = 0.7) = 1, we find 
the best numerical value for g{v4/v2 = 0.65) such that 
the corresponding flows collapse as a function of Lcff- We 
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Figure 10: As a function of temperature, for V4/V2 = —10 
. Top : Specific heat per site Cv/N. The dashed line repre- 
sents a linear interpolation of the regular part of the speci- 
fic heat. Inset : Maximum of Cv/N as a function of system 
size. Middle : Binder cumulant and inverse stiffness crossings. 
Insets : Binder cumulant derivative (at T = 0.672) and deri- 
vative L.dK~^ /dT (at T — 0.6715) as a function of system 
size (log-log scale). Bottom : Columnar susceptibility x per 
site. Inset : Maximum of x as a function of system size (log- 
log scale). All lines in insets denote power-law fits of critical 
exponents (see text for details). 



then proceed successively with the next two consecutive 
values of W4/U2 and find the best factor for g{vi/v2 = 0.6). 
In this manner, we go through all the parameter space, 
trying to collapse the curves two by two and finding the 
corresponding g{vn/v2)-, down to W4/W2 — —0.6. If such an 
action is possible, and if the estimated function g varies 
monotonically as a function of v^jv^^ then all the critical 
points in the interval [—0.6,0.7] should refer to the same 
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Figure 11: Specific heat per site C„/A^ at v^jvi = —1.0. The 
dashed line represents a linear interpolation of the regular part 
of the specific heat. Inset : Maximum of Cv/N as a function 
of system size. 



critical behavior. 

We have searched for such a scaling function g and 
we have arrived to the unambiguous conclusion that 
it is not possible to perform a global collapse of the 
fiows. In particular, the slow divergence of the flow for 
va/v2 — —0.02 is not compatible with the strong growth 
shown by the flows around W4/W2 = 0.6. This is because 
the flows between these two intervals can hardly be col- 
lapsed with their neighbors (meaning that for any pair 
of consecutive flows in the interval < W4/U2 < 0.4, 
we could not find a rescaling factor such that the two 
flows are superimposed). On the contrary, we have suc- 
ceeded in performing two local collapses : one for the 
region —0.6 < Vi/v2 < 0.02, and the other for the region 
0.4 < U4/W2 < 0.7 (see Fig. 12 bottom). This is a clear in- 
dication of the presence of two different critical behaviors 
in the phase diagram. Because the collapse for positive 
W4/W2 > is strongly diverging, we naturally associate 
it with a first order transition. The collapse in the ne- 
gative range of Vi/v2 describes most probably a second 
order regime. One can see that this collapse has actually 
not converged, as it should for a continuous transition. 
This is partly due to our original choice of the constant A 
which made our operational critical temperature slightly 
above the real T^- A small amount of non-zero stiffness 
remains present in the system and perturbs the conver- 
gence towards the fixed point behavior. 

To further confirm our conclusion on the flowgram, we 
have performed the following check. Suppose that the 
points Vi/v2 = —0.02 and W4/U2 = 0.6 refer neverthe- 
less to the same critical regime. Then, there should exist 
a scaling function connecting the flows at 114 = —0.02 
and W4 = 0.6. It is in fact possible to connect very 
roughly the two flows by renormalizing the length L for 
Vi/v2 ~ 0.6 by a factor g(0.6) = 8 while leaving the flow 
for Vi/v2 — —0.02 unchanged (see Fig. 13 left)). If this 
rescaling is really a physical renormalization transforma- 
tion, then the properties of the system at Vi/v2 = —0.02 
and L = 160 should correspond to those at V4,/v2 = 0.6 
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Figure 13; Left ; Flows at 114/^2 = —0.02 and 114/^2 = 0.6 at 
the rescaled length Left- Top right : Maximum of the specific 
heat for 114 = as a function of the system size. Data are taken 
from Ref. 16 except for the value at L = 180. The solid line is 
the power-law fit, giving rise to the estimate a/v = 1.15(15) 
compatible with the result of Ref. 16 without the knowledge 
of the L=180 point. Bottom right : Probability of energy per 
site at W4 = and L = 180. 
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Figure 12: Flowgram (top) and performed collapses (bot- 
tom) for the extended dimer model. See text for details of the 
flowgram procedure. 



and L ^ 32. For this latest value, we know in particular 
that the transition is first order since the energy histo- 
gram displays a double peak (see Fig. 5). If the critical 
point W4 = is indeed in the same regime as the critical 
point at V4/V2 = 0.6, we should expect to see a double- 
peak also for samples of sizes L > 160 at V4 = 0. We have 
therefore carried a simulation at 1)4 = with L — 180 
and have then measured the histogram of energy at the 
maximum of the specific heat. The histogram displays a 
unique and well-defined peak (see Fig. 13 bottom right), 
which confirms that the transition at W4 = is conti- 
nuous. Moreover, the value of the specific heat maximum 
is perfectly compatible with the exponents measured pre- 
viously at W4 = (see Fig. 13 top right), confirming the 
validity of the exponent a/i' ~ 1.11(15) found in Ref. 16. 
To conclude on this part, not only can we drive out the 
possibility that the regimes v/^jvi < and V4/V2 > are 
connected but we can also state that the tricritical point 
is necessarily located at a value w|/w2 > 0. 



VI. DISCUSSION 

In this study, we presented an extended version of 
the classical dimer model with plaquettes and cubic in- 
teractions. Our simulations indicate that depending on 
the sign of the interactions, the nature of the critical 
phase transition between the Coulomb and columnar 



phases changes. Compared to the pure plaquette mo- 
del {V4/V2 = 0), a cubic interaction which reinforces the 
alinement of dimers (yi/v2 > 0) leads to a first order 
transition, identified by a double peak distribution of the 
energy and a strong divergence in thermodynamic quan- 
tities. Whether an infinitesimal positive cubic coupling 
is enough to alter the nature of the transition is not cer- 
tain, as we lack any theoretical argument to support this, 
but simulations tend to indicate that this change should 
occur very close to W4 = 0. On the other side, when both 
interactions compete {V4/V2 < 0), finite size scaling ana- 
lysis of thermodynamic quantities shows no sign of any 
discontinuity at the transition. When the cubic interac- 
tion is largely dominant {vi/v2 = —10), the critical ex- 
ponents deviate significantly from those measured in the 
pure plaquette model, with a ^ 0.2, v '^ 0.6 and rj '^ 0.2 
in the first case and a ^ 0.5, v ~ 0.5 and 77 ~ in the 
latter. In the interval — 1 < V4^/v2 < 0, we find exponents 
in between these two cases, probably due to a cross-over 
effect. To discharge the possibility of a weak first order 
transition on the frustrated side, we carried out a flow- 
gram analysis in the vicinity of 1)4 = 0. The fiowgram 
clearly demonstrates the presence of two groups of flows, 
one corresponding to positive values of V4/V2 and the 
other to negative values. This rules out the possibility of 
a connection between the two parts of the phase diagram, 
and thus reveal the presence of a tricritical point^**. Fi- 
nally, we also detected the presence of a new crystalline 
phase at low temperature, deep in the frustrated regime. 
This phase is characterized by a degeneracy growing with 
the system size. 

In a parallel work, Papanikolaou and Betouras re- 
cently studied another extension of the dimer model 
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that slightly differs from ours but which leads to similar 
conclusions^^ . To perturb the pure plaquette model, these 
authors introduce further neighbor interactions between 
dimers which preserve the cubic symmetry. They find 
that when the extra couplings favor the columnar ali- 
nement, the transition becomes first order. On the op- 
posite, a frustrating coupling maintains the continuity 
of the transition. The critical exponents they measure 
are in agreement with ours in the limit of large frustra- 
tion. Our study has the advantage of using much larger 
samples (the maximum system size used in Ref. 22 is 
L = 32), giving more confidence on the values of ex- 
ponents. Furthermore, the use of the flowgram analysis 
provides an explicit evidence for tricriticality. In another 
related work, Chen et al}^ find that by favoring one par- 
ticular subset of the 6 columnar orderings at low tempe- 
rature, the transition between the Coulomb and colum- 
nar phases becomes either first order or in the id XY 
universality class. 

While not entirely interpretive, it is tempting to com- 
pare the exponents on the frustrated side of the dimer 
model to other models which also display unconventional 
phase transitions. In Ref. 18, a lattice gauge theory 
representing a coarse-grained version of the dimer model 
was simulated by Monte Carlo. There, the exponents 
were found to be very close to the 3d XY universality 
class and thus differ from those presented here. Another 



source of information is given by the study of the 
NCCP^ field theory. There exists several indications 
that the exponent rj in the NCCP^ theory should 
be large and positivc^^. For instance, the transition 
observed in the ring-exchange models studied in Ref. 26 
provide an estimate of 77 ^ 0.2 — 0.35 (depending on 
which correlation function is looked at) and v ~ 0.68(1). 
These values are not far from those we observed in the 
strong frustration regime. On the other hand, direct 
simulations of lattice versions of the NCCP^ all exhibit 
a first order transition^". This possibility is ruled out in 
the dimer model, thanks to our results obtained with 
the fiowgram analysis. 
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